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Abstract 

The invaded cluster algorithm is used to study the XY model in two and three dimensions up 
to sizes 2000 2 and 120 3 respectively. A soft spin 0(2) model, in the same universality class as 
the three-dimensional XY model, is also studied. The static critical properties of the model and 
the dynamical properties of the algorithm are reported. The results are K c = 0.45412(2) for the 
three-dimensional XY model and n = 0.037(2) for the three-dimensional XY universality class. 
For the two-dimensional XY model the results are K c = 1.120(1) and n = 0.251(5). The invaded 
cluster algorithm does not show any critical slowing for the magnetization or critical temperature 
estimator for the two-dimensional or three-dimensional XY models. 
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I. INTRODUCTION 



In this paper we present a Monte Carlo study of the XY model using the invaded cluster 
(IC) algorithm. The purpose of this study is both to obtain high precision results for the 
XY model and to test the efficacy of the invaded cluster method for phase transitions with 
continuous symmetry breaking. The XY model is in the 0(2) universality class and in two 
dimensions the transition is of the Kosterlitz-Thouless type ]T). Theoretical and compu- 
tational studies of the 0(2) critical point include high temperature series expansions [[|, 
renormalizat ion-group calculations M and Monte Carlo (MC) simulations ||, ^], || [7|. The 
A transition in 4 He is also in the 0(2) universality class and the specific heat exponent a for 
this system has been measured to very high precision & M. Ill 



Recent Monte Carlo studies of the XY transition use versions of the Wolff algorithm flT 



because near the critical point they are much more efficient than local algorithms. The Wolff 
algorithm is an example of a cluster algorithm of the kind first introduced by Swendsen and 



Wang fL2| , [L3| for Ising-Potts models. The central idea of cluster algorithms is to identify 
clusters of sites by a bond percolation process correlated to the spin configuration. The 
spins of each cluster are then independently flipped. Cluster algorithms can be extremely 
efficient when the percolation process that defines the clusters has a percolation threshold 
that coincides with the phase transition of the spins. This situation holds for the original 
Swendsen- Wang algorithm applied to Ising-Potts models and was later shown to hold for a 
variety of other spin systems with discrete symmetries [T5|, [TBI. Wolff |]1TJ showed how 
to extend cluster algorithms to spin models with 0(n)-symmetry by an Ising embedding 
method. 



Invaded cluster (IC) algorithms |L7|, [Tg, |19[ are cluster algorithms with the property that 
they find and simulate the critical point automatically without a priori knowledge of the 
critical temperature. The critical temperature is a direct output of the algorithm and the 
magnetic exponents can be obtained from finite-size scaling of the cluster size distribution. 
In addition to providing the critical temperature and magnetic exponent without having to 
invoke methods such as histogram re-weighting, IC algorithms appear to have less critical 
slowing than corresponding Swendsen- Wang or Wolff algorithms ||19|| . IC algorithms can, 
in principle, be constructed whenever a conventional cluster algorithm is available with the 
property that the bond percolation process has a percolation threshold that coincides with 
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the phase transition. IC algorithms have thus far been applied to systems with discrete sym- 
metry breaking including Ising-Potts |T7, TB, [21]] models, Widom-Rowlinson models |22 



and the fully frustrated Ising model In this paper we show how Wolff's embedding 

scheme can be used to construct an efficient IC algorithm for simulating the 0(2) critical 
point in three dimensions and the Kosterlitz-Thouless point in two dimensions. 

One of our objectives is to obtain a high precision value of the r] exponent for the 0(2) 
universality class. Because rj is itself very small, it has proved difficult to measure it with 
much precision. One of the difficulties is the presence of corrections to scaling. Recently, 
Hasenbusch and Torok || [7| proposed a method to minimize this difficulty by considering 
a soft spin 0(2) model with a parameter controlling the variance in the length of the spin 
vectors. By adjusting this parameter they minimize corrections to scaling and improve their 
estimate of rj. Below we apply a modified version of the IC algorithm to the soft spin model. 



In the next Section we describe the IC algorithm for the XY model. In Sec. |T| we give 
results for the critical temperature, magnetic exponents and dynamic properties of the XY 
model in two and three dimensions. In Sec. |IV| we describe the algorithm and present results 
for the soft spin model. Conclusions are presented in Sec. [V] 



II. INVADED CLUSTER ALGORITHM FOR THE XY MODEL 

The algorithms for the XY model used in this paper are obtained by combining the 
invaded cluster method and Wolff's embedding scheme for continuous spin models. The XY 
model is defined by the Hamiltonian: 

(3H = -K Si ■ s, (1) 

<i,j> 

where Sj is a two dimensional unit vector and the summation is over all nearest neighbor 
bonds on the lattice, here either the square or cubic lattice with periodic boundary condi- 
tions. 

We begin by describing a version of Wolff's algorithm for the XY model. On each Monte 
Carlo step we choose a two-dimensional unit vector r. For each bond (i,j) of the lattice, 
the bond is called satisfied if both spins lie on the same side of the line perpendicular to the 
unit vector, that is, 

(s t -r)(s r r)>0. (2) 
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Satisfied bonds are then occupied with probability 

P(4 sj) = 1 - exp[-2#(r • s,)(r • sj)}. (3) 

One way to implement the occupation of bonds with this probability is to independently 
assign random numbers Uij, uniformly chosen from the interval [0, 1) to every bond of the 
lattice and then to occupy the satisfied bonds if mj < P(si,Sj). An equivalent approach, 
which provides a useful link to the IC methodology, is to define Kij from u^j by substituting 
j j for P(si, Sj) and Kij for K in Eq. (|^) and then solving for K^j, 



u 



K 



'■J 



log(l-u iij )/{2{f-s i ){f-s j )). (4) 



Satisfied bonds with k^j < K are occupied. The occupied bonds define a set of connected 
clusters. Single sites with no occupied bonds are considered clusters so that every lattice 
site is uniquely a member of some cluster. 

Once the clusters are identified each cluster is flipped with probability 1/2. A cluster is 
flipped by reflecting every spin in the cluster through the line perpendicular to f, 

^ — > R(f)Si (5) 

where 

R(f)Si = Si - 2(si ■ r)r. (6) 

Flipping clusters with probability one half yields a new spin configurations and completes 
one MC step. It is straightforward to show that the algorithm is ergodic and satisfies detailed 
balance. 



Wolff's original paper [IT] introduces two innovations. The first is to grow and flip 
only a single cluster in each Monte Carlo step. The second is the generalization of cluster 
methods to spin models with continuous symmetries. Only the generalization to continuous 
symmetries is used here and combined with the invaded cluster methodology. In fact, Wolff's 
single cluster method is not compatible with the invaded cluster methodology. When we 
speak of the Wolff algorithm in the following we mean the algorithm described above for 
which all clusters are defined and flipped in a single MC step. 

Cluster algorithms can generally be viewed as a sequence of bond moves and spin moves. 
During the bond move a configuration of occupied bonds and the associated set of clusters 
are generated from the spin configuration. This is done by occupying satisfied bonds with a 



probability that depends on the simulation temperature. During the spin move a new spin 
configuration is obtained by randomly flipping clusters. The bond configurations can be 
viewed as a correlated percolation model. For a cluster algorithm to be efficient, the corre- 
lated percolation model should have a percolation transition that coincides with the critical 
point of the spin model. If this holds then clusters of all sizes are flipped during a single MC 
step and changes to the spin configuration occur on all length scales. For Ising-Potts mod- 
els the correlated percolation model associated with the Swendsen-Wang algorithm is the 
Fortuin-Kastelyn random cluster model and the equivalence of the percolation threshold and 



the critical point is well-understood [|3|, P4j . The percolation properties of the bonds defined 



by the Wolff algorithm for the XY model have recently been investigated p5| . The conclu- 
sion is the same as for the Ising-Potts models: the critical point for the three-dimensional 
XY model and Kosterlitz-Thouless point for the two-dimensional XY model coincide with 
the percolation threshold for the bonds defined by the Wolff algorithm. 

The invaded cluster methodology relies on the equivalence of the phase transition in 
the spin model and the percolation transition of the occupied bonds to find and simulate 
the phase transition without prior knowledge of the critical coupling. Like other cluster 
algorithms, a full MC step consists of a bond move followed by a spin move. The spin 
move is the same as for standard cluster algorithms but the bond moves differs. In the IC 
bond move, satisfied bonds are occupied one at a time in random order until a signature of 
percolation is first observed. The set of occupied bonds obtained in this way defines bond 
clusters and these are flipped in the usual way. The signature of percolation is incorporated 
in a stopping condition that is tested after each new bond is occupied. In this paper we use 
a topological stopping condition, which requires that at least one cluster wraps around the 
lattice in at least one direction. 

For the most part, the IC algorithm closely parallels the Wolff algorithm. First a unit 
vector is randomly chosen and satisfied bonds are determined with respect to this unit 
vector as described in the paragraph including Eq. (^|). Then uniform random numbers, Uij 
are assigned to each bond and non-uniform random numbers, obtained from them 

according to Eq. (|j). The next step of the IC algorithm differs from the Wolff algorithm. 
Satisfied bonds are occupied one at a time according to the ordering defined by Kij with the 
satisfied bond having the smallest value of k occupied first. After each bond is occupied, the 
set of clusters is updated and the stopping condition is checked. If no cluster wraps around 
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the lattice the satisfied bond with the next largest k is occupied but if some cluster wraps 
around the lattice in some direction the bond move is completed and no further bonds are 
occupied. The largest value of k chosen during the bond move is called k. The spin move is 
identical to the Wolff algorithm: with probability 1/2, each cluster is reflected through the 
line perpendicular to r according to Eqs. (|) and (|5]). 

The IC algorithm simulates the critical point and k is an estimator of the critical coupling. 
To see why this is the case, consider a large system with a spin configuration that is typical of 
the critical point. The correlated percolation threshold for the occupied bonds is related to 
critical coupling according to Eq. (§). That is, if satisfied bonds are occupied with probability 



P c (s i} Sj) = 1— exp[— 2K c (f-Si)(f-Sj)], the occupied bonds will just percolate p5| . Occupying 
bonds with probability P c (si, Sj) is the same as occupying all satisfied bonds one at a time 
in ascending order of k and stopping at the largest k such that k <K C . The IC algorithm 
works in the same way except that satisfied bonds are added until a cluster wraps around 
the system. For a large system, this event will occur with k nearly equal to K c . Thus, if the 
system is at criticality, a single step of the IC algorithm is almost identical to a single step of 
the Wolff algorithm at K c and k pa K c . So long as the fluctuations in k become small as the 
system size increase, the IC algorithm and the Wolff algorithm at K c will sample the same 
state in the thermodynamic limit. However, in finite volume, the invaded cluster algorithm 
defines an ensemble that is different than the canonical ensemble and has different finite-size 
scaling properties. 

So far we have assumed the spin system is already at the critical point. Suppose that 
the spin configuration is typical of the low temperature phase, then there will be more 
satisfied bonds than at criticality and a smaller fraction will be needed to form a spanning 
cluster so that k < K c . The invaded cluster MC step thus corresponds to a Wolff MC step 
at temperature T = J/k > T c so that the system is warmed. A similar argument shows 
that if the system starts in the high temperature phase, it is cooled by the IC algorithm 
so that there is negative feedback mechanism that forces the system to the critical point 
independent of the starting configuration. A detailed discussion of these arguments in the 
context of Ising- Potts models is given in Ref . [[Tj|] . 

For the 2D XY model, the IC algorithm described above does not perform well. The 
problem is similar to a problem that occurs in the 2D Ising model. Specifically, at criticality, 
the satisfied bonds themselves just percolate. This means that in a significant fraction of 
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spin configurations, spanning is not possible. It also means that the distribution of the 
temperature estimator is broad because of spin configurations for which spanning is just 
barely possible. The problem is even worse for the the XY model. Since the unit vector 
r is randomized from one step to the next, spanning on one Monte Carlo step does not 
guarantee spanning on the next Monte Carlo step as it does in the Ising case. As a result, 
a small fraction of the time, no spanning cluster can be found for the chosen unit vector r. 
Furthermore, the distribution of the temperature estimator is very broad because some spin 
configurations and unit vectors just barely allow spanning yielding a very small value of the 
temperature estimator. 

These problems can be alleviated by taking advantage of a second independent set of 
satisfied bonds. Let's refer to the bonds that are satisfied with respect to the definition of 
Eq. (0) as red satisfied bonds. Let b be a unit vector perpendicular to r and define bond 
as blue satisfied if 

(s l -b)(s r b)>0. (7) 

To increase the probability that we obtain a spanning cluster, we can occupy both red and 
blue satisfied bonds. Two values of k, one with respect to r and the other with respect to 
b, are assigned to each bond using Eq. (J|) and a single value of u. The entire set of ft's, for 
both red and blue satisfied bonds, is ordered. Bond are either red or blue occupied in the 
order prescribed by the k's and sets of red and blue clusters are identified and updated after 
each new bond is occupied. The first cluster to span, either red or blue, stops the process of 
occupying bonds. During the spin move, both sets of clusters are flipped with probability 
1/2 according to Eqs. and (|G|) with r replaced by b for the blue clusters. One way of 
looking at this algorithm is that it utilizes two independent embeddings of the Ising model in 
the XY model, one embedding relative to r and the other relative to b. The two embedding 
method was always able to find a spanning cluster and it produces a less broad distribution 
for the temperature estimator. Nonetheless, the distribution of k is still very broad and 
we find better statistics by averaging 1/k. All results reported for the two-dimensional XY 
model use the two-embedding method. 



7 



III. RESULTS 



The algorithm was implemented on systems with maximum linear size L = 120 for three 
dimensions and L = 2000 for two dimensions. Each run consisted of 160000 MC steps 
for three dimensions and 10000 MC steps for two dimensions . We collected statistics for 
the estimator of the critical coupling, k and the size of the spanning cluster, M. For the 
two-dimensional model we obtained the estimator of critical coupling as inverse of the 
estimator of the critical temperature. The average critical coupling and its error bars were 
obtained using the blocking method [13] with 100 blocks of 100 data points each for two 
dimensions and 160 blocks of 1000 data points each for three dimensions. The average size 



of the spanning cluster and its error bars was obtained using the bootstrap method |L3 
The reported error bars are one standard deviation. The simulation time was 10 -5 CPU 
seconds per spin per MC sweep on a 450MHz Pentium III machine running Linux. 



A. Critical temperature of the three-dimensional XY model 

Figure |1] shows the average value of k for the three-dimensional XY model versus the 
inverse of the linear size of the system 1/L. We fit the data for L > 10 to a function of the 
form: 

(«(L)> = - K \ (8) 

where K c , a and p are parameters of the fit. The results of the fit are K c = 0.45412(2), 
a = -0.64(2) and p = 1.211(9) with x 2 = 6.23, x 2 /d.o.f. = 0.69 and the confidence level 
is Q = 0.72. The value for the critical coupling K c compares well with values used in the 
literature as shown in Table [VT] . Note that p ^ 1/u « 1.5 as would be expected from naive 
finite size scaling arguments. As is the case for Ising- Potts systems, the IC ensemble for the 
XY model does not have the same finite size scaling properties as the canonical ensemble. 

The validity of the IC method is justified by the fact that the width of the distribution 
of k decreases as L increases. Figure shows the standard deviation a k of k as a function 



of 1/L and suggests that 



lim a k = (9) 

L^oo 



and therefore we expect a sharp distribution of k for L = oo. The solid curve in Fig. |2] is 
the result of a fit to the data for L > 50 to the functional form an = a + bL~ q yielding 
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a = 0.0002(10), b = 0.43(3) and q = 0.66(3). 



B. Kosterlitz-Thouless temperature 

Figure^ shows (1/k(L)) _1 as a function of 1/ log(L) for the two-dimensional XY model. 
We fit the data for L > 160 to the function, 

1 + a(log(L)) 2 

where and a are parameters of the fit. This functional form was motivated by combining 
the usual finite size scaling assumption that £ = L with the Kosterlitz-Thouless expression 
for the critical behavior of the correlation length £ f26 |, 



UK) ~ e u (11) 

where t is the reduced temperature and v — 1/2. For this fit we found K c = 1.120(1) and 
a = 2.49(4) with X 2 = 3.9, x 2 /d.o.f. = 0.65 and Q = 0.68. The choice of L > 80 was a 
compromise between keeping the error in K c small and the confidence level Q large. The 
result for K c is in good agreement with some of the recent results from the literature as 
shown in Table |V11| . Although the fit is good, the result for K c should be viewed with 
caution because it is based on finite size scaling assumptions that do not necessarily hold 
for the IC ensemble. 

Figure |3] shows the standard deviation er^ of 1/k as a function of 1/ log 2 (L). Unlike the 
situation in three-dimensions, it is not clear that aya vanishes as L — ► oo. If the finite size 
scaling behavior of the IC ensemble is of the "essential singularity" type observed in the 
canonical ensemble then a reasonable hypothesis is that ayx ~ (logL) -9 . A naive extrap- 
olation suggests aya approaches a finite value near 0.11 but a slowly decreasing function 
cannot be ruled out. In the previous section we discussed reasons for the broad distribution 
of the temperature estimator. Simulations with L ^> 2000 would be needed to determine 
whether or not ai/% vanishes as L — > oo. 

C. Magnetic Exponents 

The average mass M of the spanning cluster is proportional to the magnetization of the 
system so that the magnetic exponents can be obtained from the fractal dimension of the 



spanning cluster, defined by M ~ L . The critical exponent r/ is related to D via, 

r] = 2 + d-2D. (12) 

Figure |5| shows a log- log plot of M vs. L for the three-dimensional XY model. A linear fit 
yields r/(3D) = 0.037(2) with x 2 = 2.1, X 2 /d.o.f. = 0.35 and the confidence level Q = 0.9. 
The smallest value of L included in the fit is L min = 50. Figure || shows values of r] obtained 
from fits with smallest linear size L min . The value of rj has an upward trend up to L min = 40 
and then is constant for larger values of L m , n . The value of L m i n for the reported results was 
chosen such that the statistical error is minimal for a reasonably large Q and the value of rj 
is well into the region of constant values. The plot on Fig. |] also shows an upward trend for 
the two last values of L m i n . Table |V1] shows some recent results for the three-dimensional 
XY model. Our result for 77 agrees with recently published values. 

Figure [7] shows a log-log plot of M vs. L for the two-dimensional XY model. A linear fit 
to the data yields f](2D) = 0.251(5) which is in good agreement with the theoretical value 
for the Kosterlitz-Thouless transition t]kt = 0.25. The smallest value of L included in the 
fit is L min = 480 and x 2 — 2.84, x 2 /d-o.f = 0.94 and Q = 0.41. The value of L min was 
chosen such that the statistical error is minimal for a reasonably large Q. 



D. Dynamics of the invaded cluster algorithm 



The autocorrelation function for a time dependent variable A(t) is defined as: 



r m <(A(0)-<A>)(A(t)-<A>) > 

TA{t) ~ < (A(0))- < A >) 2 > (13) 

The integrated autocorrelation time is defined by 

r A = -+ lim £r A (t). (14) 

z t=l 

The integrated autocorrelation time for A is interpreted as the time needed to obtain sta- 
tistically independent measurements of A. As a result, the statistical error in measuring A 
in a MC simulation is proportional to (r^/iV) 1 / 2 where N is the number of MC steps. 

We measured the autocorrelation functions and the corresponding integrated autocorrela- 
tion times for the magnetization M and the critical coupling estimator k. When calculating 
the integrated autocorrelation time it is necessary to choose a finite cut-off for w in the sum 
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in Eq. (0). We used w = 100, a value much longer than the integrated autocorrelation 
times obtained below. Figure |8] shows the autocorrelation function versus time for the in- 
verse critical temperature of the three-dimensional XY model. Figure [5] shows the same 
plot for the two-dimensional case. Note the negative overshoot of the correlation functions 
for time of one MC step. This is due to the presence of a negative feedback mechanism 
in the IC algorithm. Consecutive temperature estimators for the spin configurations are 
anti-correlated. 

Table shows the integrated autocorrelation times for the mass of spanning cluster tm 
and the inverse temperature tk c as a function of the linear size of the system L in the three- 



dimensional XY model. Table IV shows the same data for the two-dimensional case. These 



data show that the autocorrelation times for M and k are independent of the the system 



size. This behavior was also observed for IC dynamics for Ising- Potts models [17, 19, 27 



and gives the impression that the IC algorithm does not have any critical slowing down. 



However, as shown by Moriarty, Machta and Chayes[19] for Ising- Potts models, observables 
defined on scales intermediate between the lattice spacing and system size have relaxation 
times that diverge with system size so that the dynamic exponent for the algorithm is greater 
than zero. Nonetheless, the small values of r k and tm mean that relatively few MC steps 
are needed to obtain good statistics for k and M. 

IV. THE SOFT SPIN 0(2) MODEL 

The values of r\ for the three-dimensional XY universality class vary significantly in the 



literature. Table |VI] shows some recent results for r\. Possible reason this discrepancy in the 
literature is the presence of finite size scaling correction in the fit of M vs. L. Hasenbusch and 
Torok || and Campostrini at. al. |7j studied a 4 soft spin model defined by the Hamiltonian: 

H/T = -J/T + £ <^ 2 + A ]>> 4 2 - l) 2 (15) 

<i,j> i i 

where T is the temperature, J interaction constant and A a "softness" parameter. The 
vectors (ft can have any value of the modulus \(f>\. The XY model is obtained as a special 
case for A = oo. This model is in the same universality class as the three-dimensional XY 
model. Hasenbusch and Torok showed that finite size corrections in the canonical ensemble 
are minimized near A = 2.0. 
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We implemented the IC algorithm for the soft spin model. Following Refs. |28| and || 



we used the IC algorithm to update the orientation of the vectors <fi and the Metropolis 
algorithm to update the modulus. The IC part of the combined algorithm is described in 
Sec. 2. The only difference is that the modulus of the spins does not necessarily have the 
value of |0| = 1. Since the IC algorithm performs only reflections, the modulus of the 
spins remains the same after the IC sweep. After every IC step, the algorithm performs an 
update of both the modulus and the orientation of the spins using the Metropolis algorithm 
as follows. New values of the two spin components are proposed: 

<t>x' = <j>x~ 2(Px - 0.5) 

<t>y = 4>y ~ 2(Py ~ °-5) 

where p x and p y are random numbers from the interval [0, 1). The update is accepted with 
the probability 

P = min[l, exp(H/T - H' /T)\ 

where T is obtained from the output of the previous IC step. The temperature will therefore 
be different for every Metropolis update. However, as is the case for the three-dimensional 
XY model, these temperature fluctuations tend toward zero as L —>■ oo. In the limit of very 
large system size the Metropolis subroutine updates the system at a fixed temperature. The 
Metropolis subroutine preserves the negative feedback mechanism of the IC algorithm. If the 
temperature of the system is low, the Metropolis algorithm will tend to make the modulus of 
the spins larger. With larger spin modulus the IC algorithm will give a higher temperature 
as an outcome. The same holds if we start with high temperature. The Metropolis algorithm 
will make the modulus of the spins smaller and consequently the IC algorithm will lower the 
temperature. 

In our simulations of the soft spin model we used A = 2.0. The maximum linear size 
of the system was L = 120 and each run consisted of 10000 MC steps. Figure [10] shows 
the average critical coupling («) vs. 1/L. A fit to the function (k) = K c /(1 + aL~ p ) yields 
K C (X = 2.0) = 0.5100(1), a = -0.9(4) and p = 1.3(1). The smallest size included in the fit 
is L = 30 and x 2 = 10.0, \ 2 /d.o.f. = 1.4 and Q = 0.18. This value of K c is in agreement 
with the value reported by Hasenbusch and Torok K c = 0.5099049(6). Figure |TI] shows 
a log-log plot of M vs. L. The resulting value for t], r](\ = 2.0) = 0.042(8) (with \ 2 = 5-8, 
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X 2 /d.o.f. = 0.98 and Q = 0.43) is in agreement with their result rj = 0.0381(2). A recent 
work by Campostrini et. yielded rj = 0.0380(4). 

V. CONCLUSIONS 

In this paper we have introduced and applied an invaded cluster algorithm for the two 
and three-dimensional XY models and a related three-dimensional soft spin model. This 
work extends the range of validity of the invaded cluster method to include continuous spins 
and systems where the phase transition is of the Kosterlitz-Thouless type. Our results for 
critical temperatures and magnetic exponents are in reasonable agreement with recent values 
in the literature. The invaded cluster algorithm is very efficient for XY systems, showing no 
critical slowing for estimators of the critical temperature and the magnetization and yielding 
highly accurate values of the magnetic exponent with relatively little computational effort. 
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FIG. 1: Critical coupling (k) vs. 1/L for the three-dimensional XY model. The solid line is a fit 
to the data as described in the text. 
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FIG. 2: Standard deviation of the critical coupling vs. 1/L for three-dimensional XY model. 
The solid line is a fit to the data as described in the text. 
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L 




R 


M 


10 


0.1079(2) 


0.4730(2) 


236.2(2) 


20 


0.0621(1) 


0.46195(7) 


1334(1) 


30 


0.04658(8) 


0.45890(5) 


3662(3) 


40 


0.03810(7) 


0.45751(3) 


7480(4) 


50 


0.03281(6) 


0.45668(2) 


13038(10) 


60 


0.02904(5) 


0.45620(2) 


20474(18) 


70 


0.02635(4) 


0.45584(2) 


30052(25) 


80 


0.02409(4) 


0.45558(1) 


41820(36) 


90 


0.02240(4) 


0.45538(1) 


56047(46) 


100 


0.02080(4) 


0.45521(1) 


72824(64) 


110 


0.01957(3) 


0.45511(1) 


92197(80) 


120 


0.01854(3) 


0.455006(9) 


114374(96) 



TABLE I: Numerical data for the three-dimensional XY model. 



L 


Urp 


k = l/T 


M 


10 


0.44(5) 


0.788(2) 


55.5(1) 


20 


0.32(5) 


0.891(1) 


200.6(5) 


40 


0.25(5) 


0.953(1) 


717(2) 


80 


0.22(5) 


0.9936(7) 


2631(6) 


160 


0.19(5) 


1.0212(8) 


9543(24) 


240 


0.19(5) 


1.0350(8) 


20151(50) 


320 


0.18(5) 


1.0410(8) 


34779(86) 


480 


0.17(5) 


1.0519(9) 


74083(186) 


640 


0.16(5) 


1.0571(7) 


126798(311) 


800 


0.17(5) 


1.0605(7) 


192100(490) 


1000 


0.16(5) 


1.0643(6) 


293519(756) 


2000 


0.15(5) 


1.0742(7) 


1073907(2800) 



TABLE II: Numerical data for the two-dimensional XY" model. 
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FIG. 3: Critical coupling l/(k x ) vs. l/ln(L) for two-dimensional XY model. The solid line is a 
fit to the data as described in the text. 
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FIG. 4: Standard deviation of the critical coupling (Ti/% vs. l/ln 2 (L) for two-dimensional XY 
model. 
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FIG. 5: Mass of the spanning cluster InM vs. InL for the three-dimensional XY model. The solid 
line is a linear fit. 
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FIG. 6: The critical exponent r\ for the three-dimensional XY model vs. minimum system size 
Lmin included in the fit of the data. The maximum system size in the fit is L max = 120. 
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FIG. 7: Mass of the spanning cluster InM vs. InL for two-dimensional XY model. The solid line 
is a linear fit. 



L 




k 


M 


10 


0.151(1) 


0.562(1) 


229(1) 


20 


0.0906(6) 


0.5486(9) 


1288(4) 


30 


0.0652(5) 


0.5160(7) 


3540(10) 


40 


0.0498(3) 


0.5139(5) 


7270(20) 


50 


0.0437(3) 


0.5151(4) 


12691(40) 


60 


0.0382(3) 


0.5123(4) 


19988(70) 


70 


0.0377(2) 


0.5120(4) 


29404(100) 


80 


0.0320(2) 


0.5116(3) 


40846(140) 


90 


0.0320(2) 


0.5115(3) 


54384(200) 


100 


0.0300(2) 


0.5113(4) 


70784(240) 


110 


0.0282(2) 


0.5111(3) 


89484(300) 


120 


0.0266(2) 


0.5112(3) 


111851(400) 



TABLE III: Numerical data for the three-dimensional soft spin model with A = 2.0. 
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FIG. 8: The autocorrelation function of the critical coupling r« vs. Monte Carlo time for the 
three-dimensional XY model. The system size is L=100. 
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FIG. 9: The autocorrelation function for the Kosterlitz-Thouless temperature estimator T 
Monte Carlo time for the two-dimensional XY model. The system size is L=1000. 
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FIG. 10: Critical coupling (k) vs. 1/L for the three-dimensional soft spin model. The solid line is 
a fit to the data as described in the text. 
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FIG. 11: Mass of the spanning cluster InM vs. InL for three-dimensional soft spin model. The 
solid line is a linear fit. 



23 



L 




T M 


10 


0.195 


0.50 


20 


0.129 


0.56 


30 


0.045 


0.52 


40 


0.052 


0.52 


50 


0.045 


0.51 


60 


0.028 


0.50 


70 


0.025 


0.50 


80 


0.021 


0.49 


90 


0.028 


0.50 


100 


0.029 


0.58 


110 


0.010 


0.52 


120 


0.032 


0.61 



TABLE IV: Autocorrelation times for critical coupling r« and magnetization tm for the three- 
dimensional XY model. 
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L 






10 


0.273 


0.73 


20 


0.144 


0.95 


40 


0.107 


0.57 


80 


0.077 


0.71 


160 


0.095 


0.88 


240 


0.061 


1.15 


320 


0.073 


0.80 


480 


0.097 


0.99 


640 


0.071 


1.07 


800 


0.082 


0.82 


1000 


0.060 


0.60 


2000 


0.073 


0.64 



TABLE V: Autocorrelation times for critical coupling r« and magnetization tm for the two- 
dimensional XY model. 
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Ref. method K c {XY ) rj 



This work (IF) MC 0.45412(2) 0.037(2 

This work (soft spin) MC - 0.042(8 

H| MC 0.45408(8) 0.036(14 

§ MC+HT - 0.0380(4 

H MC 0.454165(4) 0.042(2 

§ MC - 0.0381(2 

§ MC 0.45420(2) 0.024(6 

§ FT - 0.038(6 

|] HT 0.45419(3) 0.039(7 

@ FT - 0.0349(8 

H MC - 0.035(5 



TABLE VI: A summary of recent estimates of the critical coupling for the three-dimensional XY 
model and the exponent rj for the three-dimensional 0(2) universality class. 



Ref. method 



This work MC 1.120(1) 0.251(5) 

MC 1.118 0.238(4) 

MC 1.113(6) 

MC 1.1199(1) 0.233(3) 

MC 1.106(5) 

MC 1.1209(1) 



TABLE VII: A summary of recent estimates of the critical coupling and the exponent r\ for two- 
dimensional XY model on a simple cubic lattice. According to the Kosterlitz-Thouless theory 
rj = 1/4. 
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